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We discuss an instability in the leapfrog integration algorithm, widely used in current Hybrid Monte Carlo 
(HMC) simulations of lattice QCD. We demonstrate the instability in the simple harmonic oscillator (SHO) 
system where it is manifest. We demonstrate the instability in HMC simulations of lattice QCD with dynamical 
Wilson-Clover fermions and discuss implications for future simulations of lattice QCD. 



1. INTRODUCTION 



2. SHO RESULT AND HYPOTHESIS 



The subject of instabilities in molecular dy- 
namics (MD) integrators is not new. Studies of 
the phenomenon for the case of lattice QCD have 
been reported in In the lattices used 

were small (4"* sites). The onset of instabilities 
was investigated as a function of the integration 
step size (St) for heavy and light quark masses. 
For both masses studied the instabilities set in 
at some critical step size. The critical value of 
St was smaller for the lighter mass than for the 
heavier one. However, in both cases St was still so 
large that the energy change {SH) from the step 
size induced errors along an MD trajectory was 
sufficiently great to make simulations impractical 
in the unstable region of parameter space anyway. 

Our study using larger lattices and quark 
masses that are relatively light (by current stan- 
dards) encountered instabilities at step sizes that 
are small enough to be relevant for large scale 
simulations. Our motivation to summarise these 
results is to highlight the impact these instabil- 
ities may have on future simulations using MD 
integrators as components, in particular HMC. 

This article is organised as follows. In section ^ 
we demonstrate the instability in the simple har- 
monic oscillator (SHO) system. At this point we 
will restate the hypothesis of ^ to explain the 
onset of the instability. We present our main re- 
sults in section || and show that they are com- 
pletely consistent with this hypothesis. Finally 
in section H we will draw our conclusions. 



In the case of the SHO system a leapfrog up- 
date can be written as: 
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where q and p are the coordinates and conjugate 
momenta, t is MD time and lo is the angular fre- 
quency. 

The eigenvalues of the matrix in (pT) are 
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When luSt < 2 the inverse cosine in (||) is real, 
the eigenvalues are complex, and the phase space 
trajectories describe ellipses. However, when 
ujSt > 2, the inverse cosine becomes imaginary, 
and one eigenvalue is greater than one. At this 
point the phase space trajectories become hyper- 
bolic, and the energy change along a trajectory 
grows exponentially. 

The hypothesis of is that the short distance 
modes of an asymptotically free quantum field the- 
ory behave like a collection of loosely coupled SHO 
modes. It is then anticipated that the MD inte- 
gration of lattice QCD will go unstable when the 
highest frequency mode goes unstable. 

We note that in the SHO example uj plays the 
role of the MD force F. Hence we surmise that for 
QCD the instability will occur when the combi- 
nation FSt reaches a critical value. Furthermore, 
we expect that the fermionic contribution to the 
force will increase in some manner with the in- 
verse quark mass. This suggests that as the quark 
masses are decreased, the FSt term will reach its 
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critical value for smaller values of 6t. This indeed 
is exactly what is reported in |^ and 1^ . 

3. NUMERICAL RESULTS 

Our numerical computations were carried out 
using 10 lattice configurations from a UKQCD 
lattice computation ||3|. The lattices were picked 
evenly spaced from a larger ensemble. The lattice 
volume was 16'^ x32 sites, and the physical param- 
eters for the ensemble were /3 — 5.2 for the gauge 
coupling, K — 0.1355 for the hopping parameter 
and csw = 2.0171 for the 0(a)-improvement co- 
efficient 1^,^. These parameters correspond to 
a regime where the ratio of the pseudoscalar to 
vector masses is ^ ?s 0.6 HJ^. 

We evolved the configurations for unit length 
MD trajectories using a variety of step sizes at a 
variety of n values. Along each trajectory we com- 
puted SH as well as the 2-norms and oo norms of 
the gauge and fermionic contributions to the MD 
force (||Fg||2,||F^||2, \\Fg\\oo, and ||F/|U_) aver- 
aged over all the time steps along the trajectory. 

In figure |l| we show the results of fitting the 
fermionic force norms (averaged over configura- 
tions) to the fit ansatz: 



Fit Ansatz: |[ F || = C ( am )", am = ('4){1/k - 1/kJ 
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in order to check that the forces increase in some 
inverse manner with the quark mass. To calculate 
am we used the formula 
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where Kcrit is the critical value of k correspond- 
ing to the massless limit of the pion. In our fits 
we have left Hcrit as a free parameter, however 
to determine am for the axes of figure we used 
the value Kent = 0.13633, determined from spec- 
troscopy on the ensemble (||. 

One can see from figure l| that the values of 
a are negative for both the 2-norm and the oo— 
norm indicating that the forces do indeed grow 
inversely with am. The fits for Kcrit are consistent 
with the spectroscopic determinations indicating 
that the forces will diverge as am 0. 

In figure ^, we plot the variation of the force 
norms and 6H against the step size St at a fixed 
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Figure 1. II-F/II2 and ||i^/||oo vs am. 

value of K = 0.1355 for a single configuration. It 
can be seen that at a value of about St = 0.0105, 
the fermionic forces start increasing and that SH 
increases exponentially indicating the onset of the 
instability. The gauge contributions to the force 
norms seem to show no change in behaviour. 
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Figure 2. 2-norms (bottom), cx)- norms (middle) 
and corresponding values of SH (top) vs St. 

The 00-norms grow at a larger rate than the 
2-norm. As the oo~norm is the largest abso- 
lute value contribution to the fermionic force, it 
plays a similar role to the highest frequency mode 
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5x=0.01 and 8T-0.012, single configuration 
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Figures. \\Ff\\^ (bottom), \ \Fg\\^ (middle), SH 
(top) vs K, for two values of the step size: St — 
0.01 (circles) and St — 0.012 (squares). 

of a corresponding set of loosely coupled SHOs 
whereas the 2-norm of the force can be likened to 
the behaviour of the corresponding bulk modes. 
Hence this result is consistent with the earlier hy- 
pothesis, that the highest frequency modes drive 
the system unstable. 

In figure |^ we show the behaviour of ||-F/||oo, 
ll^gllcx) and SH against the hopping parameter k 
for St = 0.01 and St = 0.012. When St = 0.01 
the system is completely stable for all values of k 
(for this single configuration); however, for St — 
0.012 the fcrmionic force norm rises sharply for 
K = 0.1355 accompanied by a large increase in SH 
signalling the onset of the instability. This again 
is consistent with our earlier hypothesis that in- 
stabilities set in when the FSt term reaches a 
critical value. 

4. CONCLUSIONS AND DISCUSSION 

We have exhibited the onset of an instability 
in the leapfrog MD integration scheme, used in 
HMC simulations and inexact simulation algo- 
rithms. The instability is manifest in the leapfrog 
scheme for a simple SHO system. We have also 
presented numerical evidence demonstrating that 
the instability is present also in lattice QCD sys- 
tems. We have shown that our data is consis- 
tent with the hypothesis that the instability oc- 
curs when the FSt term in the leapfrog equations 



reaches some critical value, and that the fermionic 
contribution to the force increases with decreas- 
ing quark mass. We have shown that as quark 
masses become lighter the instability sets in at 
smaller values of St. Furthermore, for simula- 
tions with light dynamical quarks the instability 
sets in at values of St that are small enough to 
be relevant to large scale numerical simulations. 

We anticipate that the major stumbling block 
for simulations with lighter quark masses will not 
come from "exceptional" configurations, where 
the value of k may become supercritical for some 
configurations, but rather from the onset of the 
instability which we expect to set in before k can 
become supercritical. 

The instability problem is exponentially bad, 
but can be controlled easily by reducing St. We 
are concerned for simulations with inexact algo- 
rithms, where the instability may be hard to de- 
tect as there is no SH to monitor, and full finite 
step size extrapolations are seldom made. 

Finally we wish to mention that neither the 
use of higher order integration schemes nor the 
use of double precision arithmetic are expected 
to alleviate the instability problem 
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